A discrete time-dependent method for metastable atoms in intense fields 
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The full-dimensional time-dependent Schrodinger equation for the electronic dynamics of single-electron 
systems in intense external fields is solved directly using a discrete method. Our approach combines the finite- 
difference and Lagrange mesh methods. The method is applied to calculate the quasienergies and ionization 
probabilities of atomic and molecular systems in intense static and dynamic electric fields. The gauge invariance 
and accuracy of the method is established. Applications to multiphoton ionization of positronium and hydrogen 
atoms and molecules are presented. At very high intensity above saturation threshold, we extend the method 
using a scaling technique to estimate the quasienergies of metastable states of the hydrogen molecular ion. The 
results are in good agreement with recent experiments. 



I. INTRODUCTION 

Matter exposed to intense laser fields has attracted exten- 
sive research in the past two decades. Due to the competi- 
tion between external forces induced by the laser field and 
Coulomb interactions, which bind the system together, the na- 
ture of the metastable states that arise is strongly controlled by 
the nature of the external field and is not simply a feature of 
the atomic or molecular system. In such circumstances the 
internal and external fields are on an equal footing and not 
separable. 

The study of these effects in few-electron atoms have be- 
come accessible recently both experimentally and theoreti- 
cally because of advances in experimental techniques and 
the availability of supercomputers. Compared with atoms, 
molecules in intense laser fields are much more complicated 
not only because of the multi-center nature of the problem 
but also due to the additional vibrational and rotational de- 
grees of freedom associated with the nuclei 1 1, 2]. In poly- 
atomic systems, the energy transfer to collective motion leads 
to extremely high-energy secondary photons, ions and elec- 
trons 1 3, 4]. 

Standard perturbation theory is no longer applicable when 
the forces induced by the applied laser field are comparable 
to the binding forces of the system. If the external field is 
periodic then the system is metastable and has a well-defined 
quasienergy spectrum. Suppose the unperturbed ground-state 
has a real energy Ei, then the effect of the external field is 
to produce an energy shift A and decay width rate T/h, so 
that the quasienergy has the usual form E — Ei + A — iT /2. 
This describes a steady-state decay and can be treated very ef- 
fectively with time-independent methods such as the Floquet 
method. However, at the very highest intensities the transient 
aperiodic effects of short duration pulses of dynamic fields 
mean that the shift and rate are poorly defined. Under these 
conditions a time-dependent approach is essential. In this 
paper we describe such an approach applied to the regimes 
where the quasienegy is well-defined and to time-dependent 
problems where this is not the case. We find the method works 
efficiently and accurately for both cases, and thus is well- 
suited to the study of metastable states under all conditions. 



Recent advances in experimental technology in the utilisation 
of high intensity lasers have led to the creation and study of 
short-lived atomic and molecular states. In this paper we ap- 
ply the method to problems of this type, though the method 
is of more general applicability. With very short pulses, not 
only the laser envelope but even the phase of the field can be 
important |5]. Therefore, one has to directly integrate the full- 
dimensional Schrodinger equation in order to describe accu- 
rately the underlying physical mechanisms k^ii|8J and obtain 
data relevant to experiments. 

Although one-dimensional models have been routinely 
used in describing atoms and molecules in strong laser 
fields 0], ab initio full-dimensional quantum-mechanical 
calculations are very important for exactly calculable few- 
electron systems. As discussed in the numerical results 
for excitation, dissociation and ionization in simplified mod- 
els are strongly sensitive to the parameters chosen. Full- 
dimensionality calculations without approximation are really 
necessary to establish accurate data as benchmarks to assess 
the quality or regime of applicability of other calculations and 
to produce results which are comparable with experimental 
observations. 

We have recently developed a very accurate and efficient 
numerical method to study metastable states in high inten- 
sity fields. This was applied to the fragmentation of the hy- 
drogen molecular ion in an intense laser field beyond the 
Born-Oppenheimer approximation Q 03, fTHl . In this paper, 
we make a detailed investigation of this method applied to 
metastable atomic and molecular systems. We discuss the 
variational characteristics of this method and confirm its ac- 
curacy, convergence and gauge invariance. Results are com- 
pared to other theoretical estimates of the metastable states, 
and we compare with experimental results for the ionization 
of the hydrogen molecular ion at high intensities. 



II. DISCRETE TIME-DEPENDENT SCHRODINGER 
EQUATION 

Our goal is the direct solution of the time-dependent 
Schrodinger equation of an arbitrary one-electron systems in 
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strong external electric fields 



where 
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where the electronic potential energy is the non- separable sin- 
gular function 
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and A{r, t) is the (external) vector potential; V(r, t) is the 
(internal) scalar potential. An example of such a problem is 
a one-electron molecular system in an intense laser field; the 
internal field arising from the positively-charged nuclei. In 
order to simphfy the model and to understand the electronic 
dynamics in isolation, we consider the hydrogen molecular 
ion with the nuclei fixed in space. The nuclei, labelled 1 and 
2, are a fixed distance R apart, and have charges Z\ and Z2. 
The origin of the coordinate r is located at the internuclear 
midpoint, with the electronic Hamiltonian takes the form 



A' ZiZ2 

and the molecule-laser interaction term is 
V^\z,t)=zE{t). 
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Consider a short optical pulse that is approximately 
monochromatic, within the bandwidth limit, with a well- 
defined peak intensity. This can be simulated by choosing an 
electric field of the following form 
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with n = \r + and = \r — ^R\, and A the vec- 
tor potential. We include the constant internuclear poten- 
tial for reasons of convention. Removing the quadratic term 
{e^/2me)A^ by gauge transformation, and making the dipole 
approximation, the Hamiltonian in the Coulomb gauge can be 
written as 
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while in the length gauge: E{t) = — dA/dt 
1 
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Monochromatic light with linear polarization parallel to the 
internuclear axis implies a cylindrical symmetry about this 
axis. Associated with this symmetry is a good quantum num- 
ber, A, proportional to the projection of angular momentum 
along the axis. Thus the electron position can be completely 
described by the radial, p, and axial, z, coordinates with re- 
spect to an origin taken at the midpoint between the nuclei. 
The time-dependent Schrodinger equation reduces to a 2h-1 di- 
mensional partial differential equation. Hence in atomic units 
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E{t) = Eof{t) cos (jLt, 
where the pulse envelope, f{t), is given by 
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(10) 

and in which Eq is the peak electric field, the pulse ramp time 
is Tl and the pulse duration T2, with associated bandwidth 
Au) = 1/t2. Since the peak field, Eq, is related to the cycle- 
average intensity, (7), by the relation (7) = ^ccoEq then the 

conversion formula is _Eo ~ 5.338 x 10^^^ (7) if the intensity 
(7) is in W cm~^, and the field strength in atomic units. The 
corresponding ponderomotive energy is Up = £'o/(4a;|). 
Consider the molecule initially in the ground state X ^5]+. 
It is convenient to change the dependent variable to remove 
the first-derivative in p as follows 



Hp,z,t) = (27rp)i/V(p,z,t), 
so that the time-dependent equation is 

i^^(j){p,z,t) = [T^ + Tp + V^{p,z,R) 



(11) 



where 
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This 2+1 dimensional equation, given in Eq. ( I12t . can be 
discretized on an Np x x Nf space-time grid. We la- 
bel the Np radial nodes by, {pi, p2, . . . pi, . . . pNp}, while the 
Nz axial grid points are denoted by, {zi, Z2, ■ ■ ■ zj, . . . zjv^ }. 
The evolution progresses through the sequence of times 
{ti,t2, . . . tfc, . . . tNt}- Thus the wavefunction can be written 
as the array (t){p,z,t) (j>{pj, Zi,tk) = 4>ijk- The method 
of discretization of the Hamilton divides the axial and radial 
coordinates into subspaces. Two distinct but complementary 
grid methods are used for the subspaces. The radial subspace 
is discretized over a semi-infinite range using a small number 
Np of unevenly spaced points that are the nodes of global in- 
terpolating functions that enter the Lagrange mesh technique. 
This leads to a small dense matrix for the Hamiltonian in the 
p-subspace. On the other hand the axial coordinate subspace 
is represented by a large number of equally-spaced points that 
are the mesh points of a finite-difference scheme. The asso- 
ciated subspace Hamiltonian matrix is large but sparse. Our 
approach is tailored to the requirements of accuracy and com- 
putational efficiency. 



where these mesh functions have properties both of Lagrange 
interpolation functions and exact discrete orthogonality, that 
can be summarised as follows 
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The matrix element of an operator or function Q{x, ^), in 
this basis is given by 
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where the error, e, in the Gaussian quadrature depends on Np, 
a and the form of Q. Expanding the wave function in this 
basis 



No 



A. Discretization of tlie Hamiltonian in the radial subspace 

For wavefunction character in the p-coordinate we have 
found that a Lagrange mesh can provide a more ef- 

ficient discretization scheme compared with that provided by 
finite-difference formulae. This type of grid can be chosen 
to accommodate short-range singularities or long-range be- 
haviour, and can be scaled in length or number of grid points 
to improve accuracy with a modest number of points. We 
rescale the radial variable as follows, p — hx where h is some 
arbitrary scaling factor, and < a; < +oo. Then consider the 
set of functions 

= ( v( ^' ^.y ^^^'e-^^'Li-Hx), (15) 
\T{a + n+ 1) J 

where (x) are the generaUzed Laguerre polynomials 

(16) 
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n\ dx 

These functions form an orthonormal set on the domain 

< X < oo 



iPm{x) (pn{x) dx = S„ 
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Then, for any given value of a, one can construct an iVp-point 
grid based upon Gauss quadrature rules. Choosing the grid 
points (xi, X2, ■ ■ ■ , XNp) the Np solutions of l''^\x) — 0, 
the quadrature weights corresponding to these pivots are given 
by the Christoffel numbers A^, where 



One can define the set of differentiable functions 
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In the radial subspace the kinetic energy is represented by the 
dense matrix [12] 
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S^i = (xiXi)^^'^ ^ Xi,^{xk - Xi) ^{xk - xi) ^ (25) 

k^i,l 

Then the Hamiltonian in the radial subspace has the form 

{H,)u = Su{T,) + {Tp)u+SuV,^{p„z,R) 

+ SaV^i{z,t). (26) 

The last grid point, the largest root of L^'' (p/h) ~ 0, defines 
Pmax, the radius of the cylindrical box. 



B. Discretization of the Hamiltonian in the axial subspace 

The axial coordinate grid is chosen to be a set of N^ equally 
spaced points which cover the range — Zmax < ^ < ^max 
with a separation Az — 2zniax/{Nz — 1), so that zj = 
— z,nax + (j — l)Az. We choose the method of finite dif- 
ferences to treat this coordinate in order to make effective use 
of parallel processing and for appropriate treatment of wave- 
function dependence on z |13]. The sparsity of the matrix 
and confinement of communication to that between nearest 
neighbours is ideal for efficient calculation. For example, the 
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kinetic energy can be evaluated by the five-point finite central 
difference formula 



1 



zjjm 
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with error proportional to (Az)* and resulting in a sparse 
(pentidiagnonal) matrix. The momentum operator arises in 
the velocity gauge perturbation and is given by 



V 



ml 



-iA{pi,t) 
12(Az) 
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j+2,m + oOj+l.m 



— 8(5j_i,„i + 5j-2,m), (28) 

with error proportional to (Az)^. 

Therefore the augmented [Np x Nz) x {Np x N^) Hamil- 
tonian matrix takes the form 



il,jm 



= 5il {Tz)jm + {Tp)il Sj 



(29) 



in the length gauge and 



il,jm 



5 a {Tz)jm + {Tp)ii 5 



jm 
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Kir 
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in the velocity gauge. 



The orthonormal set is formed using Gram-Schmidt orthogo- 
nalization. Letting h^: denote the [ua + 1) x [ua + 1) upper- 
Hessenberg matrix formed by the coefficients (hfc)y we ob- 
tain the matrix equation 



hfc = QiHfcQfc, 



(36) 



where is a matrix formed from the Ua column vectors 
[qQ, qj^, . . . , q„^] and so is the Krylov subspace Hamil- 
tonian which is calculated simultaneously with Q^,. V\k = 
QjjhfeQj. can be used as a replacement to in a wide variety 
of applications. In particular the time-evolution operator can 
be written as 



Qfee 



(37) 



Now h/j is typically a tridiagonal matrix and so e^*'^''^* can 
be computed inexpensively. Thus \}{tk + Ai, tk) is in effect 
a unitary propagator correct to order (At)"" . 

III. ACCURACY OF THE DISCRETE SOLUTIONS 

The results we obtain from this method depend on several 
parameters, namely Np, hp, Az, Zmax, na. At, and the choice 
of gauge used to describe the laser-molecule interaction term. 
In this section we detail the choices for these parameters re- 
quired to obtain accurate and fully converged solutions. 

A. Spatial grid parameters 



C. Discretization and propagation in time 

The wavefunction is discretized on a grid so that it makes 
up a vector of Np x Nz components. At time tk the (i, j)th 
element of the vector can be defined as 

{Mk)i,] ^ (t>{hxi,Zj,tk), (31) 
while the (zZ, jm)th element of the matrix is defined 

(Hfe).,,,>„ = [H{tk)]a,,m , (32) 
so that the time evolution is described by the equation 

H(t)v(t) = iv(t). (33) 

Suppose that time is divided so that tk+i = t^ + At then the 
solution can be propagated using the unitary evolution matrix 

Vfc+i = U(tfc + At, tk)vk ~ expi-iHk At)vfe. (34) 

Evaluation of this exponential is carried out using a Krylov 
subspace decomposition 0, Using the Arnoldi al- 

gorithm fl4ll we construct an orthonormal set of vectors, 
[qp, q;^, q2, . . . , q„^] , which span the Krylov subspace 

i^„„(Hfc, Vfe) = span {vfc, HfeVfe, H^v^, . . . , H^"Vfe} . 

(35) 



In an intense field problem, the dynamics of the system can 
be very sensitive to the initial state. Therefore, the accuracy 
of the ground state wavefunction and its corresponding energy 
play a crucial role. Several important grid parameters are open 
to choice, and to determine these we proceed as follows. We 
first choose the values of Az and Np, then apply a variational 
method by adjusting the value of < hp < 2 until we get 
an accurate ground-state energy. Importantly, once this pro- 
cedure is completed for a specific internuclear distance R, the 
same scaling factor works well for all R. In this hybrid finite- 
difference/discrete-variable method we noted that there is a 
delicate relation among the values of Az, Np and hp. For a 
fixed Np value, if we double the value Az, we have to roughly 
double the value of hp in order to maintain accurate energies. 
Since the kinetic energy terms are homogeneous in these coor- 
dinates this maintains the balance between these terms. In the 
following calculations, we take the range of the z coordinate 
as [-Zmax, Zmax], whcrc Zniax = 300.9 a.u. and Az = 0.1 
a.u. (6019 points in all). The z-subspace is shared across the 
processor array, in our case we used 13 processors for this 
task. The p-subspace is spanned fully on each processor and 
we take Np = 30 with hp = 0.5185 corresponding to the limit 
Pmax ~ 55 a.u. 

To obtain the ground state eigenvectors it is convenient and 
efficient to use an iterative Lanczos method fs^. InTableHI we 
list our ground state energies for A = and A = 1 at different 
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TABLE I: Calculated energies for the lowest states of gerade sym- 
metry with A = and A = 1 for different 7?, compared with the 
exact values from Sharp's tabulation 115,1 . Values of R and energy 
are in atomic units. 



R 



A = state 
Exact Present 



A = 1 state 
Exact Present 



1.0 -0.451785 -0.451783 0.525893 0.525872 

2.0 -0.602635 -0.602636 0.071229 0.071216 

4.0 -0.546085 -0.546088 -0.100825 -0.100830 

6.0 -0.511968 -0.511972 -0.130325 -0.130327 

8.0 -0.502570 -0.502574 -0.134511 -0.134512 

10.0 -0.500580 -0.500582 -0.132716 -0.132716 

12.0 -0.500167 -0.500172 -0.129950 -0.129948 

16.0 -0.500035 -0.500040 -0.126253 -0.126243 

20.0 -0.500015 -0.500018 -0.125084 -0.125072 



internuclear distances and compare them with the exact values 
from Sharp's tabulation. Here, we have taken Az = 0.1 a.u. 
and Np = 30 in both calculation. However, the scaling param- 
eter hp has been adjusted to 0.5185 for A = and 0.146 for 
A = 1, these values are used for all R in the table. The results 
are in excellent agreement compared with the exact results. 

The excited state spectrum supported by the grid can be 
found using the spectrum of the autocorrelation function ^ 



Cit) 



dVr(r,t)0(r,O), 



(38) 



where (/)(r,0) is an arbitrary trial function, and 4'{r,t) is 
the function evolved in the field-free Hamiltonian (t>{r,t) = 
exp{—iHot/h)(j>{r,0). The natural frequencies (eigenener- 
gies) appear as peaks in the spectral density 



P(^,T) 



(39) 



The resolution improves with longer times, T. The trial func- 
tion can be chosen to find the states of given symmetry. We 
have compared the energies of the electronic states at an equi- 
librium separation R = 2.0 a.u. with the exact results given 
by Sharp 1 15]; this comparison was previously performed by 
Dundas |8]. Our current results are more accurate and better 
resolved than those given in |8]. This has been achieved using 
a finer axial grid Az = 0.1 a.u. as compared with Az = 0.2 
a.u. |8|, but also by adjustment of the scaling factor hp. By 
trial and error, we determined that hp 0.5185 gave the best 
estimates over the full spectrum. The scaling factor can be 
considered a variational parameter |12]. This is a useful and 
reliable method to determine the optimal grid parameters for 
subsequent dynamic calculations. 

At the very least, the dimensions of the cylindrical box, 
height 2z„iax radius pnmx, must be chosen to encompass the 
tightly-bound states of the system. However, the evolution 
through highly-excited diffuse states and low-energy contin- 
uum states is crucial to the ionization mechanism at low fre- 
quencies. The rescattering mechanism, by which slow photo- 
electrons are driven back to the core region by the laser field. 




20 30 
Time (fs) 



FIG. 1: Population loss for different time steps At. The charac- 
teristics of the laser pulse are, A = 800 nm, J = 3.2 x lO" W 
cm~^, with Ti = 5 cycles (13.4fs) and T2 = 10 cycles (26.7 fs) 

and the bond length is 7? = 5 a.u. The curves correspond to: 

At = 0.02 ; At = 0.04; At = 0.06 in atomic units. 



means that the box size should be large enough to allow the 
continuum states to evolve unfettered. For example, the free- 
electron classical amplitude of displacement and momentum 
are proportional to Eo/u}\ and Eq/ujl, respectively and so 
box dimensions of several hundred atomic units may be re- 
quired. 

The boundary between bound and free electrons was estab- 
lished in two ways. We have defined an inner region which 
satisfies either 



or 



+ {z~ R/2Y < n, 



+ (z + i?/2)2 < 



(40) 



(41) 



where R is the internuclear separation distance and rjnncr is 
taken to be 20 a.u. We also define an outer box population to 
be the total population within the entire grid. 

In order to prevent reflection of ionizing wavefunction from 
the edges of the grid wavefunction splitting which acts in the 
same manner as an absorbing potential is applied in both the 
z and p directions |8|] . 



B. Time propagation parameters 

A 12-th order Arnoldi propagator (ua = 12) was used in 
these calculations. A test of the method for numerical stability 
with respect to the time step. At using this choice of propaga- 
tor order, na, in conjunction with the spatial grid parameters 
outlined above is given in figure [2 The logarithmic scale ac- 
centuates the differences in the residual population at the end 
of the pulse. In this case we note the very good convergence 
of the method. Results are presented for the following case: 
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Time (fs) 



FIG. 2: Electron population within the inner and outer boxes for 
length (full lines) and velocity (dotted lines) gauges. The two curves 
to the left correspond to the inner box, those to right, to the outer 
box. The bond length is _R = 6 a.u., and the laser pulse parameters 
are A = 800 nm and / = 3.2 x lO" W cm~^ with n = 5 cycles 
(13.4fs) and T2 = 10 cycles (26.7fs). 

A = 800 nm and / = 3.2 x 10^'' W cm-'^, and with n = 5 cy- 
cles (13.4 fs) and T2 ~ 10 cycles (26.7fs) with R ~ 5 a.u. On 
this figure we note that only the largest time step At = 0.06 
a.u. produces a small but discernible deviations in the results. 
We find that At ~ 0.05 a.u. is sufficient for convergence in 
most cases, however the range At — 0.01 — 0.03 a.u. provides 
more reliable and accurate results. 



C. Gauge invariance 

In Fig. 12] we show the populations within the inner and 
outer boxes in both length gauge and velocity gauge. The 
bond length is i? = 6 a.u., and the laser pulse parameters are 
A = 800 nm and / = 3.2 x lO^"* W cm^^ with n ^ 5 
cycles and T2 = 10 cycles. Excellent agreement between both 
gauges is observed. This agreement also extends to longer 
wavelengths and longer bond lengths. 



IV. RESULTS 

A. Static field ionization rates 

The static field ionization rates of the hydrogen atom and 
molecular ions are known to a high degree of accuracy us- 
ing time-independent methods and provide an important test 
of our method. In our time-dependent approach the field is 
switched on over a time ti — 2 fs, with the static field main- 
tained at a constant value F for T2 = 6 fs. The rise in the field 
should be slow enough to ensure an adiabatic transition for the 
field-free ground state to the metastable state. In table |ll| we 
compare our results for the atom field ionization rates with the 



TABLE II: Static-field ionization rates F for the hydrogen atom. 
Comparison between time-independent (Floquet) calculations fT^I 
and time-dependent methods (this work). The electric field strength 
F is given in atomic units. The ionization rates are quoted in the 
format a(n) = a x 10" fs"^. 

F(a.u.) 0.1 0.08 0.06 0.05338 0.04 
Floquet 0.601(0) 0.188(0) 0.2I3(-1) 0.664(-2) 0.162(-3) 
Present 0.600(0) 0.188(0) 0.2I3(-1) 0.664(-2) 0.I63(-3) 

highly accurate time-independent results . Very good 

agreement is obtained in all cases. 

The same method can be applied to the molecular ion. 
However in this case the comparison is not so straightforward. 
Plummer and McCann 1 16] noted that at large bond lengths 
the nearly degenerate pair of Sg „ field-free states are strongly 
split by the external field. The correlated eigenstates are a pair 
of localised atomic resonances with large differences in their 
energies and their widths 71 , 72. If we apply a static field over 
a time ti — 5 fs, relatively short in comparison to that charac- 
terising the gerade-ungerade splitting, that is the hopping time 
for the electron betweeen the centers, then the electron divides 
equally between the atom sites, creating an equal mixture of 
the resonance states rather than an adiabatic transfer to one 
or other state. Since 71 ^ 72 this would imply that half of 
the population is trapped. Consider the duration of the static 
field Ts such that 71 ^ T^^, then the population decay in 
the time-dependent model is given by F = —dP/dt « ^71, 
where 71 is the width of the short-lived state. If this were the 
case, then at longer times Ts ^ 7f ^ but Ts <C 7^^ the pop- 
ulation should reach a limit of 0.5. Figure|3lshows the popu- 
lation as a function of time withri = 5 fs and T2 = 390 fs for 
F — 0.04 and R = 11. The population decays gradually to a 
limiting value of 0.5 as predicted. However when the data is 




400 



Time (fs) 

FIG. 3: Electron bound-state population as a function of time for the 
ion. Static electric field strength F — 0.04 a.u., intemuclear dis- 
tance J? = 11 a.u.. At longer times (> 200fs) the decay is inhibited 
by a trapped state created during the rise of the electric field over a 
time n = 5fs. 
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FIG. 4: Ionization rate F of Hj as a function of internuclear distance 
i? in a static electric field F = 0.04 a.u. The results show that for 
the large bond lengths, the electron splits evenly into a long-lived 
trapped state (g) and a short-lived resonance (u) : , Floquet cal- 
culation L16.1 for the u state; • present calculation x 2; O , present 
calculation. 



plotted on a logarithmic scale we notice that the decay process 
is not purely exponential. At the beginning, between 8 fs and 
13 fs, we estimate the ionization rate to be 8.75 x lO^'^ fs^^ 
and near the end 3.91 x IQ-^ fs'^ between 390 fs and 395 
fs. The ionization rate calculated by the Floquet method is 
1.76 X 10-2 fs-i for the u-state and 2.57 x 10"^ fs^^ for the 
g-state fid]. The factor of 2 difference is consistent with the 
wavefunction splitting and population trapping. Further evi- 
dence is provided in figure |3 with ionization rates as a func- 
tion of R, compared with those results obtained by Plummer 
and McCann 1 16] for the rate 71. At larger R values the split- 
ting is almost exactly one half in agreement with the adiabatic 
trapping model. However, at smaller values of R we note that 
the degeneracy of the molecular states is removed and the ion- 
ization drops rapidly as the electron is able to move adiabati- 
cally into the trapped state. 



B. Energy shift of metastable states 

The method can be applied to calculate the real part of the 
quasienergy, that is the Stark shift of the levels. In this case, 
we calculate the autocorrelation function C{t) for a trial func- 
tion evolving with the external field on. We first calculate the 
shifts for the hydrogen atom and compare with other well- 
established results. In table |lll| results are compared with 
those obtained by the complex-coordinate approach 1 18| for a 
static field DC Stark effect. Very good agreement is achieved. 
The resolution in the time-dependent method is limited by 
the bandwidth theorem: Aco w 1/Tp, where Tp is the du- 
ration of the pulse. The calculation of the AC Stark shift is 
done in the same way. We compare with the Floquet method 
applied to the hydrogen atom 1 1 8 1 for the angular frequency 



TABLE 111: Stark shifts A for the hydrogen atom ground state. The 
electric field amplitude Fq and laser angular frequencies are given in 
atomic units, ujl ~ corresponds to the static field. 



Fo 


0.10 


0.0354 







0.375 


Present results 


-0.02746 


-0.0117 


Floquet method [18] 


-0.02742 


-0.0119 



uiL ~ 0.375. This is exactly resonant with the state and 
thus, for moderate or low intensities, the AC Stark shift is to a 
good approximation half the Rabi frequency for the transition, 
that is 



A; 



1 



(42) 



where (pis and (f>2p^ are the Is and 2pz wavefunctions re- 
spectively for the hydrogen atom, Fq is the maximum of 
the electric field strength. For — 0.0354 a.u. then 
A K, —0.0132 a.u. A more precise estimate using our code 
is —0.0117 a.u., which is in good agreement with the result 
of Maquet et al 1 18] A = -0.0119 (see table|II3. Consider 
now the quasienergies of the molecular ions. In figure |5l the 
spectral density for i? = 2 a.u. is given for the laser pa- 
rameters A = 800 nm, / = 5 x 10^^ W cm-^. In this 
case, the periodicity of the Floquet spectrum is clearly vis- 
ible: i?i + A ± riLdL. The periodicity is extended and the 
peaks become sharper as the duration of the pulse increases. 
However, the resolution is limited by broadening due to the 
ionization process. The gap between any two neighbouring 
peaks is exactly one photon energy, i.e., 0.057 a.u. for the 
present case. From the data we estimate the Stark shift for 

—0.018 a.u. and for the first ex- 



the ground state is A^ 



cited state A„ 



-0.0075 a.u. The quasienergy spectrum 




-1.0 



-0.8 -0.6 -0.4 -0.2 
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FIG. 5: The quasienergies of in the presence of an intense laser 
field. The power-spectral density of the autocorrelation function is 
shown for a range of angular frequencies. In this case A = 800 nm, 
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Frequency (a.u.) 



FIG. 6: Spectral density of the correlation function showing the 
quasienergy spectrum. Comparison of the quasienergies of in the 
presence of intense laser field using different gauges. Here,A = 800 

nm, / = 5 X 10^^ W cm^^ and R = 2 a.u.. , length gauge; 

- - - -, velocity gauge. 




Laser Intensity (W/cm ) 



FIG. 7: Ionization probability of Positronium at different laser inten- 
sities for wavelength A = 780nm and pulse duration Tp = 50 fs. O 
Present Calculations; ■ Results of Madsen et al. (l^ . 



results obtained using the velocity-gauge are also shown in 
figure |6l Both spectra have the same spectral line structure 
domain shown in the figure. The method is quite useful for 
calculation of isolated resonance shifts. However when over- 
lapping resonances are present, for example at larger values 
of R the resultant spectrum is very unclear and the method 
breaks down. 



C. Ionization of positronium 

The code can be easily adapted to the ionization of positro- 
nium in an intense laser field. An investigation of this kind 
has been made by Madsen et al fioll using a time-dependent 
basis-set expansion. The scaling factor hp is adjusted to be 
0.52085 for Ps with Az = 0.1 a.u. and Np = 30 to optimize 
the ground state energy to —0.250000001 a.u.. For the sake 
of comparison 1 19], the velocity gauge is used and the pulse 
is taken as 



A{t) = Aosin^ ( ^ ) sin(t^Li), 



where Aq = Fq /lol with Fq the peak electric field strength. 
We take a pulse duration Tp = 50 fs and a wavelength A = 
780 nm. The ionization probability versus laser intensity is 
shown in figure^ Our results are in very good agreement with 
those of Madsen et al 1 1 9J over a wide range and especially for 
the intermediate intensities. 



D. 



Ionization rates for H J by intense infrared liglit 



At very high intensities the bound states will ionize ex- 
tremely quickly with a non-exponential decay. While the ion- 
ization rate or width F is mathematically well-defined, one 



cannot calculate the rate so easily from observation of popu- 
lation decay. In physical terms, in an experiment this corre- 
sponds to saturation, that is the molecules fully ionize before 
the pulse has finished. In this case it is difficult for experi- 
ments to analyse the response of the system. This requires 
extremely short pulses, with associated rapidly-varying pulse 
envelope and broad bandwidth. Under such circumstance a 
time-dependent treatment is indispensable. The fragmenta- 
tion of by an intense infrared laser (A — 790 nm) is a 
problem of current interest. For an intensity / = 3x 10^^ W 
cm~^, and R = 5 a.u. the ionization rate is roughly 5.5fs~^. 
Thus the molecule will be fully ionized within a fraction of 
an optical cycle. Mathematically, the ionization rates in this 
regime can be calculated by scaling as shown by Madsen et 
al. 1 19|. Using the arbitrary dimensionless parameters a and 
f3 we can change the scale of length and time according to 



f — af3r, i — a0^t, 



(44) 



(43) so that 



UJL- 



(45) 



The corresponding TDSE becomes 



i_a 

2^ \dz'^ dp^ p dp 



+ Ve{R,p,z) 



ip{p,z,i), 



(46) 
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TABLE IV: Test of the scaling procedure. Ground-state energies and 
the corresponding ionization rates are calculated using scaled (a — 
1.43, l3 — 1.43) and unsealed Hamiltonians. The laser parameters 
are A=800 nm and I =3.2 x lO'^* W cm~^ and a range of bond lengths 
R are considered. Energy is in atomic units, ionization rate in fs^^ 
and R in atomic units. 

Scaled model Unsealed model 



R Energy Rate Energy Rate 



4.0 -0.54600 0.128 
5.0 -0.52442 0.168 
6.0 -0.51185 0.341 
7.0 -0.50532 0.299 
8.0 -0.50232 0.289 
9.0 -0.50112 0.211 



■0.54608 0.138 

■0.52442 0.165 

■0.51199 0.351 

■0.50559 0.294 

■0.50257 0.264 

■0.50119 0.196 



with 



Zl 



Z2 



(3^~p^ + {z + R/2Y 
ak^ Z1Z2 



and 



The laser intensity and the ionization rate scale as 



and 



(47) 



(48) 



(49) 



(50) 



accordingly. 

The numerical stability scaling can be tested at intensities 
below saturation. In table II VI we compare the scaled ground- 
state energies and ionization rates obtained with those from 
a direct calculation for A = 800 nm and / = 3.2 x IQ^^ W 
cm^^ at various internuclear distances. The agreement is rea- 
sonable with relative eiTors in ionization rates below 9%. The 
values of the scaling parameters in this case are a = 1.43 and 
13 — 1.43. The interest in obtaining rates at high intensity fol- 
lows experimental work that estimated the bond length depen- 
dence of ionization rate from the ion energy spectrum li20ll2lll . 
Our results are compared with experimental spectra in fig- 
ures |8l and |9] taking a — 1.43 and (3 — 1.43. The experi- 
mental measurements collect ions over a large part of the fo- 
cal volume, and hence we adjust the normalization of the data 
to the theoretical results. It is surprising and remarkable to 
note that both spectra, experimental and theoretical, are dom- 
inated by a single large peak and that its location is repro- 
duced accurately by the theory. For R < 7 the shape of the 
ionization rate is quite well reproduced. This is surprising in 
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Internuclear distance (a.u.) 

FIG. 8: Ionization rates as a function of internuclear distance (R) 
with A = 790 nm and / = 3 X 10^^ W cm . theoretical 
calculations using Hamiltonian rescaling procedure; • , experimen- 
tal measurements |2I]. The experimental data are normalized to the 
theoretical result al R = 5 a.u.. 



view of the fact that the theoretical model is greatly simplified 
and does not include nuclear vibrations nor the averaging over 
molecular orientation and focal volume as would be required 
for a true comparison. One conclusion might be that the ex- 
tremely good agreement indicates that the ionization rate is 
a strongly-peaked function of bond length, molecular orien- 
tation and laser intensity so that the averaging process does 
not broaden these features. For internuclear distance greater 
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FIG. 9: Ionization rates as a function of internuclear distance iR) 
with A = 800 nm and / = 1.4 x 10^^ W cm"^. A, theoretical 
calculations using Hamiltonian rescaling procedure; ■, theoretical 
modelling with population depletion based on Coulomb explosion 
included, dP/dR ~ —C{V/v)P, where we take the constant C = 
0.2; O , experimental measurements |2G]. The experimental data are 
normalized to the theoretical result at J? = 5 a.u.. 
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than 7 a.u., there is no dependence on bond length, indicating 
the loss of molecular effects. An essential assumption in the 
theoretical model is that the molecules are equally populated 
at all R corresponding to the ionization process occurring as 
the molecule dissociates at steady speed. An explanation of 
the experimental shortfall in ions from large R is that if the 
molecule ionizes fully at smaller R it cannot survive to yield 
ions at large bond lengths and the ion yield rapidly drops. The 
other point in figure |8] deserving note is the overall leftward 
displacement of our theoretical curve compared with the ex- 
perimental one. There is the possibility that the experimen- 
tal calibration of the laser intensity underestimates the actual 
intensity experienced by the molecules. A very rough simu- 
lation of the molecule depletion is shown in figure|9lis which 
the theoretical ion yield is exponentially attenuated. The curve 
shown is P{R) = T{R) exp(-C /^^^ T{R)/v{R) dR) where 
v{R) is the relative velocity of the ions taken to be the clas- 
sical value for Coulomb repulsion + I/-R0 = 
The factor C is an empirical constant taken to fit to the ex- 
perimental curve in figure|9] We find C ^ 0.2 gives the best 
shape for the distribution. Taking C — 1 leads to a severe 
loss of ions at large R, almost no ions survive beyond R = 7 
in this approximation. Of course one might expect that the 
low intensity focal averaging process might raise the yield of 
ions from large R and give a more realistic picture of the pro- 
cess. This would require inclusion of the attenuation correc- 
tions discussed above. Nonetheless, it is fair to compare this 
theory and experiment for small R values, and in this respect 
the agreement is remarkably good. 



V. CONCLUSIONS 

We have made a detailed investigation of a method which is 
designed to solve the reduced-dimensionality time-dependent 
Schrodinger equation for metastable systems in intense fields. 
We have checked the reliability of the present code by exam- 
ining the convergence and the gauge dependence. Applica- 
tions to several problems have been carried out and yield good 
agreement with other available theoretical results. However, 
by direct solution of the TDSE, our method can be applied 
to both short and long pulses and to a large variety of wave- 
lengths. The provision of parallel computer architecture offers 
the opportunity to study such systems from first principles and 
in full dimensionality. 
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